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In  aquatic  ecosystems,  primary  production  (the  transformation  of  inorganic  materials  and 
light  into  living  matter  by  photosynthesis)  is  operated  mainly  by  small,  unicellular  algae 
that  float  freely  in  the  upper  layers  of  oceans  and  lakes  and  are  collectively  called  phytoplank¬ 
ton ,  see  for  an  illustration  the  phytoplankters  depicted  in  figure  1.  Since  phytoplankton 
need  light,  they  are  confined  to  the  water  layer  where  solar  radiation  can  penetrate.  This 
region  is  called  the  “euphotic  layer”  and  it  has  a  maximum  depth  of  about  one  hundred 
meters.  In  the  ocean,  this  roughly  corresponds  to  the  depth  of  the  mixed  layer;  thus,  the 
environment  where  phytoplankton  live  is  a  highly  energetic  fluid  region  characterized  by 
intense  turbulent  mixing. 


Figure  1:  Some  examples  of  phytoplankton  cells 


Aquatic  ecosystems  are  characterized  by  the  essential  role  played  by  fluid  dynamics. 
The  small  organisms  which  compose  the  plankton  are  advected  by  the  surrounding  flow 
and  must  cope  with  environmental  currents,  turbulence,  and  waves.  And  those  organisms 
which  anchor  themselves  to  the  rocks  and  to  the  bottom  must  prevent  being  wiped  away  by 
turbulent  bursts.  Fluid  dynamics  determines  the  motion  of  the  individual  organisms  and,  on 
a  larger  scale,  the  distribution  and  concentration  of  entire  populations.  Fluid  dynamics  and 
turbulence  enter  marine  and  lacustrine  Ecology  at  all  levels,  and  must  be  properly  studied 
and  understood  in  order  to  obtain  a  correct  description  of  aquatic  ecosystem  functioning. 
On  the  other  side,  the  transport  and  advection  of  a  biologically  reactive  tracer,  such  as 
the  plankton  distribution,  or  of  individual  small  living  “impurities”  (as  we  could  call  them) 
in  a  turbulent  flow  make  fluid  dynamics  even  more  intriguing,  opening  up  new  problems 
and  challenges  which  can  stimulate  the  interest  of  the  more  mathematically-oriented  types. 
In  these  three  lectures,  we  shall  touch  upon  some  of  these  issues,  with  no  attempt  at 
completeness  but  rather  barely  scratching  the  surface  of  this  topic. 
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1  Introduction 


Plankton  distributions  in  the  ocean  and  lakes  are  highly  inhomogeneous.  One  of  the  reasons 
for  this  inhomogeneity  lies  in  the  presence  of  strong  horizontal  advection,  associated  with 
mesoscale  turbulence,  coupled  with  the  spatial  inhomogeneity  of  the  nutrient  input  and  the 
nonlinear  dynamics  of  plankton  populations. 

Phytoplankton  are  ‘primary  producers’,  which  means  that  they  take  inorganic  materials, 
such  as  nitrogen  and  carbon,  and  convert  them  into  biomass  via  photosynthesis.  The  main 
limiting  factors  for  phytoplankton  growth  are  light  and  nutrient  availability;  for  this  reason, 
phytoplankton  populations  are  confined  to  the  upper  layers  of  lakes  and  oceans,  in  a  region 
known  as  the  euphotic  layer,  see  figure  1.  In  the  ocean,  nutrients  are  supplied  by  different 
mechanisms:  (1)  upwelling  of  deeper  waters,  the  most  crucial  process  as  the  majority  of 
nutrients  are  generated  by  bacterial  activity  at  depth;  (2)  direct  bacterial  regeneration  of 
nutrients  from  dead  biomass  and  excreted  organic  substances  in  the  euphotic  layer;  and 
(3)  dust  deposition  on  the  surface.  As  photosynthesis  becomes  less  efficient  when  the  solar 
radiation  is  too  intense,  and  since  most  nutrients  come  from  deeper  waters,  phytoplankton 
tend  to  live  near  the  middle  of  the  euphotic  layer. 


Phytop. 

pop 


Nutrients 


Figure  2:  The  structure  of  the  upper  ocean,  as  modeled  in  this  lecture.  Phytoplankton 
populations  are  confined  to  the  euphotic  layer  (turquoise),  where  there  is  enough  sunlight 
for  photosynthesis.  Nutrients  are  supplied  mainly  by  upwelling  of  deeper  waters  into  the 
euphotic  layer,  as  the  highest  concentration  of  nutrients  is  at  depth,  where  sea  snow  tends  to 
settle  and  bacterial  activity  is  strong.  Too  much  sunlight  is  detrimental,  and  phytoplankton 
populations  tend  to  be  found  near  the  middle  of  the  euphotic  layer. 


We  model  meso-  and  large-scale  interactions  between  fluid  dynamics  and  biology  by 
resorting  to  reaction-advection-diffusion  equations.  The  reaction  terms  represent  biologi¬ 
cal  interactions.  The  advection  terms  here  represent  horizontal  advection,  mainly  due  to 
mesoscale  circulations  and  flows.  The  diffusion  terms  are  a  parametrization  of  small-scale 
turbulence  and  of  vertical  mixing  due  to  up-  and  down-welling. 
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2  React  ion- Diffusion  Equations 

Reaction-diffusion  equations  have  the  form 

^  =  f  (P)  +  DV2P,  (1) 

where  p  is  the  concentration  field  of  some  (vector)  quantity,  and  D  is  diffusivity  (assumed 
constant  and  the  same  for  all  field  components).  The  vector  quantity  p  represents  the  con¬ 
centration  of  different  scalar  fields  undergoing  chemical  and/or  biological  interactions.  The 
vector  quantity  f  (p)  represents  the  reactions  between  the  various  components  of  p. 


In  one  spatial  dimension  and  for  a  scalar  field,  this  equation  becomes 


dp 

dt 


f(p)  +  D 


d2p 
dx 2  ’ 


where  x  is  the  spatial  coordinate  and  now  p  is  a  scalar. 


(2) 


The  simplest  biological  example  is  the  KISS  model  for  plankton  blooms  [7,  17].  This 
includes  a  term  describing  Malthusian  (exponential)  growth  and  a  diffusion  term, 


dp 

dt 


vp  +  D 


d2p 
dx 2  ’ 


(3) 


where  v  is  the  exponential  growth  rate  of  a  spatially  homogeneous  population.  Via  Fourier 
analysis,  we  can  find  the  following  expression: 


1 

p(x,t)  = —  /  p(k)  exp(ffex)  exp[(zv  —  Dk2)t]  dk,  (4) 

2vr  J_ oo 


and  so  there  is  instability  (growth)  for  k2  <  v/D  (large  spatial  scales),  while  growth  is 
restrained  at  small  scales.  This  behavior  leads  to  the  presence  of  patchiness  in  the  concen¬ 
tration  field  p  with  the  minimum  scale  2ir(D /v)1^2  ■ 


To  saturate  the  exponential  growth  of  the  homogeneous  population  at  large  times,  one 
can  include  an  extra  term  in  p2,  and  obtain  the  celebrated  Fisher  equation: 


dp 


dt 


P 


-£  =  "P  l-£  +D^, 


K 


d2P 


dx2 ' 


(5) 


where  I\  is  a  constant.  We  can  make  this  equation  non  dimensional  by  rescaling  space, 
time  and  the  concentration  field: 


P  , 

p^p  =  w.  * 


t  =  1st,  x  — ►  X  = 


v  \  V2 
DJ  ' 


88 


and  we  obtain,  omitting  the  tilde  for  ease  of  notation, 

dP  (,  v  ,  d2p 

M=p{1-p)  +  W  <6) 

where  all  quantities  are  now  non-dimensional.  The  homogeneous  version  of  this  equation, 
for  dp/dx  =  0,  is  known  as  the  logistic  equation, 

w-*1-')-  <7) 

The  Fisher  equation  admits  a  traveling  wave  solution  of  the  form  U(z)  =  p{x  —  ct ), 
which  satisfies  the  following  equation: 

?E  +  cW+Uil-U)  =  0,  (8) 

that  is, 

£  =  f =  ~cV-U(l-U)  (») 

where  (U,V)  =  ( U,dU/dz ).  This  system  of  ordinary  differential  equations  defines  a  two- 
dimensional  dynamical  system,  which  has  two  fixed  points,  namely  (U,  V )  =  (0,  0)  and 
(U,V)  =  (1,0).  Linear  stability  analysis  of  the  fixed  points  gives  the  eigenvalues  A±: 


(0,0) 

A±  =  i  -c±  (c2  -  4)1/2  , 

(10) 

(1,0)  - 

A±  =  ^  -c±  (c2  +  4)1/2  , 

Thus,  (0,0)  is  a  stable  node  if  c2  >  4,  and  a  stable  focus  if  c2  <  4;  the  fixed  point  (1,0) 
is  a  saddle.  A  stable  focus  would  result  in  negative  values  of  U  on  some  trajectories  near 
(0,  0),  and  thus  c2  >  4  is  the  only  physical  solution.  The  corresponding  phase-plane  picture 
is  shown  in  figure  3,  along  with  a  depiction  of  the  form  of  the  traveling  wave  solution  (from 
[10]). 

For  general  initial  conditions  of  the  form 

p(x,  0)  =  1,  X  <  X\\  p(x,  0)  =  0,  X  >  X2,  (11) 

where  p(x\  <  x  <  X2,0)  is  an  arbitrary  function  which  joins  1  and  0,  at  large  time  the 
solution  will  tend  to  the  traveling  wave  solution,  U(x  —  ct)  with  c  =  2,  see  figure  4. 

If  we  replace  the  second  initial  condition  as  follows 

p(x,  0)  «  Aexp(— ax),  x  — *  oo,  (12) 

see  figure  4,  then  the  system  tends  to  the  solution 

p(x,  t)  «  Aexp(— a[x  —  ct]),  x  — >  oo,  (13) 

where 

c  =  a- 1 —  if  0  <  a  <  1;  c  =  2  if  a  >  1. 
a 
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Figure  3:  Panel  (a)  shows  the  phase  plane  trajectories  for  the  dynamical  system  (9),  rep¬ 
resenting  the  traveling  wave  solution  of  the  Fisher  equation  (6),  and  panel  (b)  shows  the 
traveling  wave  solution  of  equation  (6),  for  c  >  2.  Both  illustrations  are  from  [10]. 


Figure  4:  Qualitative  description  of  the  initial  conditions  which  lead  to  a  traveling  wave 
solution,  U(x  —  ct ),  of  the  reaction-diffusion  equation  (6),  where  panel  (a)  shows  (11)  and 
panel  (b)  shows  (12).  The  linear  form  of  the  initial  conditions  for  x\  <  x  <  x 2  adopted  here 
is  for  illustration  purposes  only. 
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3  Reaction- Advection-Diffusion  Equations  and  Plankton  Dy¬ 

namics 

In  a  fluid  flow,  reactions  are  affected  not  only  by  diffusion,  but  also  by  advection.  In  such 
conditions,  the  equations  for  the  concentrations  of  reactive  tracers  become 

^  +  u-Vp  =  f(p)  +  DV2p,  (14) 

where  u  is  the  advecting  flow  field  and  D  is  the  diffusivity  (again  assumed  constant  and  the 
same  for  all  components  of  the  concentration  field). 

Often,  the  fluid  velocity  u  is  turbulent,  and  this  is  certainly  true  for  the  case  of  large 
and  mesoscale  ocean  flows.  Turbulent  advection  in  aquatic  ecosystems  occurs  on  several 
scales.  Some  of  these  are: 

Micro-scales:  Between  about  1  mm  and  a  few  meters,  the  interactions  of  3-D  turbulence, 
buoyancy  and  individual  plankton  cells  are  important,  as  will  be  discussed  in  the 
following  lecture. 

Meso-  and  submeso-scales:  Between  about  1  km  and  200km,  largely  in  the  form  of 
fronts  and  vortices,  possibly  associated  with  strong  vertical  velocities,  which  modulate 
the  plankton  distribution. 

Large  scales:  Transport  at  ocean  basin  scales  by  large-scale  gyres. 

In  the  ocean,  mesoscale  advection  is  highly  inhomogeneous,  with  very  localized  upwelling 
regions  leading  to  inhomogeneous  nutrient  fluxes;  figure  3  shows  examples  of  mesoscale 
structures  in  the  ocean  dynamical  properties  and  in  the  plankton  distribution.  In  general, 
horizontal  advection  is  associated  with  much  larger  velocities  than  vertical  advection,  and  is 
responsible  for  transporting  nutrients,  phytoplankton  and  zooplankton  as  (almost)  passive 
tracers.  Advection  timescales  in  the  upper  ocean  are  on  the  order  of  1-7  days,  which  are 
similar  to  the  phytoplankton  and  zooplankton  growth  timescales,  which  are  of  the  order  of 
1-2  days  and  1-2  weeks,  respectively.  This  means  that  scale  separation  between  the  reaction 
and  the  advection  time  scales  (which  could  allow  for  some  simplifications)  is  not  possible. 

In  keeping  with  the  tradition  of  considering  idealized  situations,  we  introduce  a  simpli¬ 
fied  description  of  mesoscale  turbulence  in  the  ocean.  Mesoscale  ocean  flows  are  dominated 
by  rotation  and  are  usually  characterized  by  stable  stratification  (ignoring  the  localized 
regions  of  intense  deep  convection  in  some  polar  areas  and  in  the  Gulf  of  Lyons).  In  such 
conditions,  one  of  the  simplest  (certainly  too  simple)  descriptions  of  mesoscale  ocean  flows 
is  provided  by  quasigeostrophic  turbulence  in  its  barotropic  version,  and,  discarding  the 
effects  of  the  free  surface  and  of  differential  rotation  (the  /3-plane  effect),  by  the  formulation 
known  as  two-dimensional  (2D)  turbulence.  Of  course,  this  flow  model  cannot  be  taken 
as  a  true  representation  of  mesoscale  ocean  flows  but  rather  as  an  idealized  picture  that 
has  some  of  the  properties  of  the  real  flow  (for  example,  the  presence  of  intense  coherent 
vortices) . 
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Figure  5:  Examples  of  mesoscale  structures  in  plankton  populations  and  in  the  ocean  tem¬ 
perature  and  velocity  fields.  Panel  (a)  shows  the  concentration  of  chlorophyll-a  around 
Tasmania  [11].  Panel  (b)  shows  a  phytoplankton  bloom  of  Nodularia  spumigena  in  the 
Baltic  sea  [12].  Panel  (c)  shows  the  intricate  mesoscale  structures  in  sea-surface  temper¬ 
ature  in  the  North  Atlantic  [13].  Panel  (d)  shows  the  velocity  field  (arrows)  and  the  sea 
surface  height  variations  (colorscale)  in  a  patch  of  the  Mediterranean  Sea,  near  the  coast  of 
France  and  Italy  [2]. 
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(15) 


The  equation  of  motion  for  2D  turbulence  is  written  as 

d(  d(  d( 

777  +  — I"  v^r  —  D  +  F, 

at  ox  ay 

where  D  is  a  (hyper)diffusion  term,  F  is  forcing,  £  is  vorticity, 

,  dv  du  _ n  , 

^  =  V' 


(16) 


and  ^  is  the  stream  function.  The  horizontal  velocities  are  given  by  u  =  —d'tjj/dy  and 
v  =  dijj/dx.  Many  theoretical  and  numerical  explorations  of  this  system  have  shown  that  a 
random  initial  vorticity  field  organizes  into  a  collection  of  coherent  vortex  structures  which 
contain  most  of  the  kinetic  energy  and  enstrophy  (squared  vorticity)  of  the  flow  (see  [9] 
for  a  pioneering  work  on  this  subject).  Figure  6  shows  an  example  of  the  vorticity  field 
produced  by  numerical  integration,  with  a  pseudo-spectral  code,  of  the  forced  and  dissi¬ 
pated  2D  turbulence  equation  (15),  see  [14]  for  further  details.  Here  we  use  hypeviscosity 
at  small  scales  and  large-scale  dissipation  to  balance  the  forcing,  obtained  by  keeping  fixed 
the  power  spectrum  at  a  given  wavenumber  kf.  After  a  transient,  the  flow  settles  into  a 
statistically  stationary  state  where  the  average  energy,  enstrophy  and  number  of  coherent 
vortices  oscillate  around  a  constant  value. 


Figure  6:  Vorticity  field,  £,  from  the  numerical  simulation  of  equation  (15)  with  forcing  and 
dissipation,  showing  the  presence  of  intense  coherent  vortices. 


4  The  role  of  localized  upwelling  regions 

Plankton  growth  depends  on  the  availability  of  nutrients,  which  have  to  be  supplied  to  the 
oceanic  mixed  layer  to  balance  biological  consumption.  In  turn,  the  nutrient  input  fluxes 
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Figure  7:  Green  regions  indicate  where  the  intense  nutrient  flux  takes  place.  The  left  panel 
refers  to  the  case  with  a  single  region  of  strong  upwelling  (RF1),  while  the  right  panel  shows 
a  snapshot  of  the  case  where  the  active  regions  are  fragmented  (RFn,  here  n=13).  The  total 
area  of  the  regions  with  strong  upwelling  is  approximately  the  same  in  the  two  cases. 


are  affected  by  different  physical  mechanisms  such  as  turbulent  mixing,  Ekman  pumping, 
convection,  transport  from  the  boundaries  etc.  Such  mechanisms  act  on  a  variety  of  scales 
providing  spatially  and  temporally  localized  fluxes  of  nutrients  to  the  biological  system. 
Thus,  it  is  useful  to  understand  the  behavior  of  a  biological  model  in  response  to  such  in¬ 
homogeneous  and  intermittent  forcings. 

The  plankton  ecosystem  model  that  is  used  here  includes  three  components,  that  repre¬ 
sent  nutrient,  N,  phytoplankton,  P,  and  zooplankton,  Z  (see  [16]  for  details).  The  dynamics 
of  the  ecosystem  is  described  by  the  reaction-advection-diffusion  equations 

~Dt  -  -  P1^TNP+ 

+m  ((1  -  +  upP  +  LizZ2)  +  dv2n 

t£=  PisHnP  -  -  f  pP  +  <17) 

‘m  =  +  »V2Z  . 

where  D/Dt  =  d/dt  +  ud/dx  +  vd/dy. 

The  terms  on  the  right  hand  side  of  the  equation  for  the  nutrient  represent  respectively 
vertical  nutrient  supply  from  deep  water,  nutrient  consumption  by  phytoplankton,  bacterial 
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regeneration  of  the  dead  organic  matter  into  nutrients  in  the  euphotic  layer,  and  diffusion. 
Bacterial  processes  are  here  assumed  to  be  faster  than  the  other  processes;  consistently, 
we  do  not  model  bacterial  dynamics  explicitly  and  assume  that  there  is  a  large  bacterial 
population  which  adjusts  instantaneously  to  the  availability  of  dead  matter  to  decompose. 
All  variables  are  measured  in  unit  of  nitrogen  concentration. 

Phytoplankton  dynamics  is  regulated  by  production,  dependending  on  available  nutri¬ 
ents  through  a  Holling  type-II  (or  Michaelis-Menten)  functional  response,  by  a  Holling  type 
III  grazing  by  zooplankton,  by  linear  mortality  and  by  diffusion  (see  for  example  [8]  for  an 
overview  of  mathematical  models  in  Ecology  and  their  jargon).  Finally,  zooplankton  grows 
when  phytoplankton  is  present  (7  is  the  assimilation  efficiency  of  the  zooplankton) ,  and  has 
a  quadratic  mortality  term  used  to  close  the  system  and  parameterize  the  effects  of  higher 
trophic  levels. 

In  the  nutrient  equation,  the  term  /i?v  is  smaller  than  one  and  represents  the  fact  that 
not  all  biological  substance  is  immediately  available  as  nutrient:  the  fraction  (1  —  /j,v)  is 
lost  by  sinking  to  deeper  waters.  Note,  also,  that  nutrient  enters  this  model  by  affecting 
the  growth  rate  of  phytoplankton.  Since  the  formulation  adopted  here  is  two-dimensional 
in  the  horizontal  and  no  vertical  structure  of  the  fields  is  allowed,  vertical  upwelling  has  to 
be  parameterized. 

Nutrient  is  brought  to  the  surface  from  deep  water  by  (isopycnal  and  diapycnal)  turbu¬ 
lent  mixing  and  upwelling.  When  these  processes  are  sufficiently  intense,  it  is  reasonable  to 
assume  that  the  surface  water  becomes  saturated  in  nutrients  (with  respect  to  the  nutrient 
content  in  deep  water)  and  further  mixing  does  not  change  the  concentration  of  available 
nutrients.  Only  when  nutrient  is  removed  (by  phytoplankton  growth  and/or  horizontal  dis¬ 
persion),  vertical  mixing  becomes  effective  again.  This  situation  can  be  represented  by  a 
relaxation  flux,  where  the  concentration  of  nutrient  at  the  surface  relaxes  to  a  value  that 
depends  on  the  nutrient  content  in  the  deep  reservoir.  The  nutrient  supply  can  then  be 
written  in  restoring  form, 

=  ~s(x,y)  (N  -  N0)  ,  (18) 

where  No  is  the  (constant)  nutrient  content  in  deep  waters  and  s  is  the  (spatially  varying) 
relaxation  rate  of  the  nutrient,  which  is  large  in  regions  of  strong  vertical  mixing  and  small 
in  regions  of  weak  vertical  mixing.  This  form  can  also  be  interpreted  as  the  finite-difference 
approximation  to  a  vertical  turbulent  advective  term  that  acts  between  two  layers  with 
nutrient  concentration  N  and  No,  and  it  is  the  standard  formulation  used  for  chemostat 
models  when  the  reservoir  has  infinite  capacity. 

Advection  is  simulated  by  the  action  of  a  barotropic  2D  turbulent  flow  -  an  approxima¬ 
tion  to  mesoscale  turbulence  in  the  ocean.  We  consider  the  case  of  statistically  stationary 
2D  turbulent  field,  forced  at  the  non-dimensional  wavenumber  k  =  40  (that  is,  at  a  length 
scale  which  is  l/40th  of  the  domain)  and  having  resolution  of  5122  grid  points.  Assuming 
that  the  forcing  scale  corresponds  to  the  typical  size  of  an  eddy,  of  about  25  km,  then  the 
domain  size  becomes  1000  km  in  dimensional  units  and  the  resolution  is  about  2  km.  The 
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turbulent  velocity  field  has  mean  eddy  turnover  time  Te  =  2.8  days. 


The  nutrient  sources  are  represented  by  a  given  number  of  randomly  distributed,  circu¬ 
lar  regions  with  fixed  radius  where  nutrient  input  is  stronger  than  average.  In  our  simplified 
view,  we  simulate  the  strong  nutrient  flux  in  these  patches,  which  we  call  “active  regions”, 
by  a  value  of  the  relaxation  rate,  s  =  sa,  which  is  100  times  larger  than  the  relaxation  rate 
in  the  rest  of  the  domain,  which  we  call  the  “passive  region”  and  where  s  =  sp  with  sp  <C  sa. 

In  all  the  numerical  simulations  discussed  here,  the  turbulent  velocity  field  is  the  same. 
What  changes  from  one  simulation  to  another  is  the  spatial  arrangement  of  the  region  where 
the  flux  is  strong  (i.e.,  of  the  “active”  region).  While  the  total  area  of  the  active  region  with 
strong  upwelling  is  kept  constant  at  the  12%  of  the  domain  area,  the  intense  upwelling  is 
respectively  confined  to  a  circular  patch  at  the  center  of  the  domain  (RF1),  or  to  a  number 
of  small  patches  randomly  distributed  in  the  domain  (RFn,  where  n  indicates  the  number 
of  individual  patches).  The  two  types  of  active  regions  are  shown  in  figure  7. 

To  characterize  the  system  behavior,  we  estimate  the  mean  primary  production,  defined 
as  PP  =  (/ 3NP/{kn  +  N)),  where  the  angular  brackets  indicate  average  over  the  whole 
domain.  Primary  production  can  be  taken  as  an  indicator  of  the  efficiency  of  the  biological 
model  to  convert  inorganic  into  organic  matter. 

The  horizontal  stirring  induced  by  the  turbulent  velocity  field  displaces  parcels  of  water 
that  are  rich  in  nutrient  and  planktonic  life  into  areas  of  limited  upwelling,  where  the  supply 
is  not  sufficient  to  sustain  the  biological  activity  at  the  level  present  in  the  parcel.  Vicev- 
ersa,  parcels  with  poor  nutrient  content  and  limited  planktonic  abundance  can  be  displaced 
into  active  regions  where  the  newly  available  nitrate  can  stimulate  a  plankton  bloom.  This 
effect  depends  on  the  intensity  of  turbulent  stirring.  The  exchange  rate  of  water  parcels 
between  active  and  inactive  regions  depends  on  the  flow  characteristics  and  on  the  amount 
of  parcels  that  are  close  to  the  boundary  areas  between  active  and  inactive  regions,  where 
there  are  strong  gradients  in  the  biogeochemical  properties  of  water.  The  more  fragmented 
are  the  upwelling  regions,  the  larger  are  the  boundary  areas.  Experiments  RF1  and  RFn  are 
designed  to  illustrate  the  effects  of  a  variation  in  the  stirring  rate  on  the  biological  activity 
when  the  size  of  the  boundary  areas  with  strong  biogeochemical  gradients  is  changed  from 
small  (RFl)  to  large  (RFn),  while  the  total  upwelling  area  is  kept  constant. 

The  enhanced  stirring  induced  by  advection  increases  the  mean  flux  from  deep  waters. 
The  enhanced  flux  originates  at  active  locations  when  a  parcel  of  water  that  has  low  nutrient 
content  is  advected  over  them.  To  see  how  this  happens,  consider  two  nearby  parcels:  one  is 
in  a  region  with  small  nutrient  upwelling  and  characterized  by  a  steady-state  nutrient  con¬ 
centration  N* ;  the  other  is  in  an  active  region  and  has  a  steady-state  nutrient  concentration 
N*.  In  this  configuration,  the  total  nitrate  flux  associated  with  these  two  parcels  of  water  is 
[sp(No  —  N* )  +  sa{N()  —  TV*)).  Suppose  now  that,  due  to  advection,  the  two  parcels  switch 
their  position:  the  parcel  with  small  nutrient  content  gets  in  a  strong  upwelling  region  and 
viceversa.  In  this  configuration,  the  vertical  flux  is  ( sp(No  —  N*)  +  sa(No  —  N*)') .  The  net 
variation  of  the  nutrient  flux  is  (sa  —  sp)(N*  —  N*),  which  is  proportional  to  sa  —  sp.  This 
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Figure  8:  On  the  left  is  the  time  evolution  of  the  average  primary  productivity.  Different 
curves  represent  simulations  with  different  number  of  fragmented  sources  (from  1  to  7). 
On  the  right  is  the  dependence  of  the  long-term  primary  productivity  (that  is,  the  value 
achieved  after  the  transient  of  panel  a)  on  the  radius  of  the  active  regions.  The  total  area 
of  the  regions  with  strong  upwelling  is  approximately  the  same  in  the  two  cases. 


term  is  positive  as  larger  relaxation  rates  are  found  in  regions  with  strong  vertical  mixing. 
The  enhanced  nutrient  flux  is  thus  due  to  the  asymmetry  in  the  relaxation  times  between 
the  active  and  inactive  regions.  Note,  also,  that  the  exchange  rate  of  water  parcels  between 
the  two  types  of  region  directly  affects  the  increased  nutrient  flux  to  the  surface,  determin¬ 
ing  larger  values  of  the  nutrient  flux  in  the  case  RFn  than  in  the  case  RF1. 

The  results  reported  here  indicate  that  primary  production  has  a  strong  dependence  on 
the  fragmentation  level  of  the  localized  forcing,  which  is  in  turn  associated  with  the  presence 
and  the  properties  of  mesoscale  structures.  These  results  imply  that  the  functioning  of 
marine  ecosystems  is  sgnihcantly  affected  by  the  flow  structure  and  by  the  distribution  of 
the  upwelling/downwelling  regions,  which  is  determined  by  the  nature  of  mesoscale  and 
submesoscale  turbulence. 

5  Advection  over  topography 

As  discussed  above,  primary  productivity  in  the  ocean  depends  on  the  availability  of  nutri¬ 
ents  in  the  euphotic  layer,  which  in  turn  depends  on  the  presence  of  upwelling  which  brings 
nutrients  from  deeper  waters. 

The  properties  of  upwelling  are  determined  by  the  structure  of  the  vertical  velocities. 
Vertical  velocities  are  greatly  affected  by  the  variations  of  topography  which  take  place  at 
steep  continental  shelves  and  near  seamounts  -  areas  where  large  abundances  of  plankton 
and  fish  are  usually  observed.  In  the  following,  we  discuss  the  impact  of  a  vortex  propagat¬ 
ing  across  a  steep  topographic  slope  on  the  dynamics  of  the  plankton  population  [19]. 

The  biological  model  is  the  same  NPZ  model  used  in  the  previous  section,  with  the  only 
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Figure  9:  A  cyclonic  vortex  on  a  sloping  bottom  topography.  Vorticity  is  on  the  left; 
divergence  on  the  right. 
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Figure  10:  Contours  of  horizontal  divergence  and  N,  P  and  Z  fields  at  t=8  days.  Biological 
fields  are  plotted  as  deviations  from  their  equilibrium  values. 
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Figure  11:  Time  series  of  N,P  and  Z  and  flow  divergence  for  two  points  at  depth  450  meters. 
Solid  line:  N,  dashed:  P,  dashed  dotted:  Z. 
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difference  that  now  the  vertical  nutrient  supply  is  given  by  the  term  <£at  =  —  (sv  +  sn)(N  — 
Nq),  where  the  rate  sv  is  proportional  to  the  vertical  velocity  of  the  flow  and  sn  represent 
a  weak  positive  nutrient  supply  in  normal  conditions  that  allows  for  sustaining  a  non-zero 
population  at  all  times.  In  the  proximity  of  a  vortex,  sv  3>  sn. 

The  value  of  sv  is  determined  by  a  Quasi-2D  dynamical  flow  model  in  the  presence 
of  topography  [4],  The  flow  is  set  in  motion  by  initially  placing  a  cyclonic  vortex  on  the 
escarpment,  see  figure  9.  The  vortex  is  affected  by  the  topography  and  it  drifts  leftward, 
according  to  a  topographic  /3-effect,  generating  a  set  of  propagating  topographic  Rossby 
waves  associated  with  a  pattern  of  vorticity,  divergence  and  vertical  velocity  with  alternat¬ 
ing  sign.  Such  motion  provides  an  oscillating  source  of  nutrients  to  the  ecological  model. 

Figure  10  shows  a  snapshot  of  the  flow  divergence  and  of  the  N,  P  and  Z  fields  after 
8  days,  when  the  topographic  wave  is  travelling  along  the  coast.  The  time  series  for  N,  P 
and  Z  and  the  flow  divergence  are  shown  in  figure  (11)  for  two  points  located  in  the  sloping 
region  at  a  depth  of  450  m  and  spaced  100  km  from  each  other.  The  divergence  time  series 
clearly  shows  an  oscillatory  pattern  with  a  period  of  about  11  days.  The  amount  of  nutri¬ 
ents  increases  whenever  there  is  positive  divergence  (that  is,  upward  vertical  velocities)  and 
decreases  if  it  is  negative  (that  is,  a  downward  velocity).  As  a  response  to  this  oscillating 
flux  of  nutrients,  oscillatory  behavior  is  observed  also  in  the  phyto-  and  zooplankton  concen¬ 
trations.  This  suggests  that  point  measurements  of  ecological  fields  over  steep  slopes  could 
reflect  the  characteristics  of  topographic  waves  rather  than  internal  ecosystem  dynamics;  a 
proper  interpretation  of  the  measurements  should  thus  take  into  account  the  properties  of 
the  Rossby  waves. 

6  Mesoscale  Turbulence  and  Coexistence 

As  a  final  example,  we  discuss  how  mesoscale  vortices  affect  species  coexistence.  To  study 
this  problem,  we  consider  a  system  with  one  limiting  nutrient,  N,  and  two  phytoplankton 
populations,  P±  and  P2  (see  [1,  15]  for  further  details).  The  system  is  assumed  to  be  in  a 
chemostat,  i.e.  there  is  a  constantly  replenishing  supply  of  N  and  a  constant  loss  of  Pj : 


DN 
~Dt 
DP i 
Dt 


-so(N-No) 


1  PiNPi 

pi  ki  +  N 


1  P2NP2 

P2  ^2  T  N 


+  DV2N, 


PiNPj 

h  +  N 


PiPi  +  DV~Pi  i  —  1,2. 


(19) 

(20) 


The  term  —sq(N  —  Nq)  is  the  nutrient  input,  described  as  a  relaxation  to  the  nutrient 
concentration  in  the  deep  reservoir,  Aro,  with  a  rate  so-  The  terms  represent  nutrient 

uptake  by  phytoplankton,  using  a  Michaelis-Menten  (or  Monod)  functional  form.  The  con¬ 
stants  pi  and  P2  are  used  to  transform  phytoplankton  biomass  into  nutrient  concentration  if 
they  are  described  by  different  units  (for  example,  nitrogen  concentration  for  nutrient  and 
carbon  or  biomass  for  phytoplankton).  The  terms  PiPt  are  natural  plankton  mortalities 
(which  include  possible  sinking  of  the  plankton  cells)  and  the  last  terms  represent  diffusion. 
Again,  the  total  time  derivative  is  D/Dt  =  d/dt  +  ud/dx  +  vd/dy.  No  direct  bacterial 
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regeneration  term  has  been  included. 


In  homogeneous  conditions  at  equilibrium,  there  is  no  co-existence  of  the  two  plankton 
populations,  and  only  the  most  favored  species  survives,  as  determined  by  the  values  of  the 
different  parameters.  In  fact,  the  so-called  Principle  of  Competitive  Exclusion  [3,  5]  states 
that  if  two  species  are  too  similar,  they  cannot  coexist  in  equilibrium:  whenever  two  species 
compete  for  the  same  resource,  the  most  favoured  will  survive  and  the  less  favoured  will 
eventually  go  locally  extinct.  By  extension,  at  equilibrium  the  number  of  species  compet¬ 
ing  on  the  same  resources  cannot  be  larger  than  the  number  of  resources.  Phytoplankton, 
however,  seem  to  escape  this  limitation,  since  a  large  number  of  species  that  compete  for 
the  same  few  resources  is  usually  observed.  This  phenomenon,  known  as  the  Paradox  of  the 
Plankton ,  was  formulated  by  Hutchinson  about  fifty  years  ago  [6].  Competition  is  avoided 
in  many  instances  by  partitioning  space  and/or  time:  the  unfavoured  species  may  be  segre¬ 
gated  in  a  spatial  environment  forbidden  to  the  stronger  species,  or  the  two  species  might 
perform  differently  at  different  times.  Mesoscale  vortices  can  be  one  of  the  causes  of  the 
spatial  segregation  of  unfavoured  and  favoured  competitors,  and  the  sheltering  effect  offered 
by  the  vortices  can  allow  unfavoured  competitors  to  survive  for  prolonged  periods  of  time. 

To  study  this  issue,  we  consider  system  (20)  with  initially  inhonogeneous  plankton  dis¬ 
tributions  and  consider  two  alternatives  for  the  flow  advection: 

1.  A  stochastic  process  representing  a  random  walk  with  memory;  the  various  ‘kicks’  are 
correlated.  This  flow  has  no  coherent  structures. 

2.  2D  turbulence,  with  coherent  vortices. 

For  the  stochastic  process,  we  use  a  Ornstein-Uhlenbeck  (Langevin)  process  [18],  defined 
as: 
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(21) 


where  TjJ  is  the  velocity  autocorrelation  time,  a  is  the  standard  deviation  of  velocities,  and 
Wx ,  Wy  are  white- noise,  Gaussian  Wiener  processes  characterized  by 


(dW)  =  0 

(dW(t)dW{t'))  =  S(t  —  t')dt, 


(22) 

(23) 


such  that  the  noise  intensities  at  different  times  are  not  correlated.  Given  these  assumptions, 
the  autocorrelation  of  the  velocity  components  is 


R(t) 


(u(t)u(t  +  t)) 

(um 


exp(— t/Tl), 


(24) 
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Figure  12:  Panels  (a)  and  (b)  show  snapshots  from  numerical  simulations  of  two  competing 
species,  in  white  and  black,  described  by  equations  (19,20),  at  the  beginning  (when  they  are 
saparate  in  two  different  regions)  and  at  a  later  time,  when  advection  has  transported  them 
across  the  flow  domain.  The  left  panels  in  both  a)  and  b)  refer  to  advection  by  the  Ornstein- 
Uhlenbeck  process  described  in  equation  (21),  and  the  right  panels  refer  to  advection  by 
2D  turbulence.  Panel  (c)  shows  the  concentration  of  the  less  favoured  population  over 
time  for  the  Ornstein-Uhlenbeck  process  (dotted  line)  and  for  2D  turbulence  (solid  line). 
The  most  favoured  population  ‘wins’  in  both  cases,  but  the  vortices  of  2D  turbulence  act 
as  transport  barriers  and  slow  down  the  contact  of  the  two  populations  and  the  rates  of 
biological  reactions. 


102 


such  that  velocities  are  correlated  on  a  timescale  T^.  The  probability  distribution  of  the 
velocity  components,  p(u)  =  p{y )  since  the  flow  is  statistically  isotropic,  is  given  by 


p(u) 


vk exp  {-h)  ■ 


(25) 


Figure  12  shows  the  comparison  between  the  plankton  concentrations  obtained  with  the 
two  types  of  advecting  flows,  simulated  in  the  same  doubly  periodic  domain  and  with  the 
same  kinetic  energy  and  decorrelation  times,  comparable  with  those  of  mesoscale  turbulence 
in  the  upper  layer  of  the  ocean.  The  most  favoured  population  (here  represented  in  white) 
‘wins’  in  both  cases,  but  it  takes  much  longer  in  the  case  of  2D  turbulence.  In  general, 
in  2D  turbulence  the  vortices  act  to  segregate  the  populations  and  allow  the  less  favoured 
population  to  survive  longer,  providing  a  shelter  area.  In  case  the  two  populations  have 
different  fitnesses  in  different  seasons,  the  presence  of  the  vortices  allows  for  the  survival  of 
the  temporarily  less  fit  population  and  for  a  long-term,  seasonally  oscillating  coexistence 
between  the  two  plankton  populations.  In  conclusion,  these  results  indicate  that  advection 
significantly  affects  biological  reaction  rates,  and  so  it  must  be  taken  into  account  when 
studying  the  dynamics  of  marine  ecosystems. 


Conclusion 

In  this  Lecture  we  saw  that  the  spatial  and  temporal  distribution  of  plankton  at  scales 
between  a  few  km  and  a  few  hundred  km  is  heavily  affected  by  the  properties  of  mesoscale 
(and  subnresoscale,  but  we  did  not  explore  it  here)  turbulence.  If  there  is  a  message  to  be 
extracted  from  these  examples,  it  is  that  the  functioning  of  aquatic  ecosystems  is  largely 
dependent  on  the  fluid  dynamics  and  that  turbulence  cannot  be  ignored. 
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